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Correcting for Bias of Kinetic Confinement Parameters Induced by Small Time Series 
Sample Sizes in Single-Molecule Trajectories Containing Measurement Noise 

Christopher P. Calderon "Q 

^ Numerica Corporation, 4^50 Hahns Peak Drive, Loveland, Colorado, 80538 

(Dated: April 16, 2013) 

Several single-molecule studies aim to reliably extract parameters characterizing molecular con- 
finement or transient kinetic trapping from experimental observations. Pioneering works from single 
particle tracking (SPT) in membrane diffusion studies [Kusumi et al., Biophysical J., 65 (1993)] ap- 
pealed to Mean Square Displacement (MSD) tools for extracting diffusivity and other parameters 
quantifying the degree of confinement. More recently, the practical utility of systematically treating 
multiple noise sources (including noise induced by random photon counts) through likelihood tech- 
niques have been more broadly realized in the SPT community. However, bias induced by finite time 
series sample sizes (unavoidable in practice) has not received great attention. Mitigating parameter 
bias induced by finite sampling is important to any scientific endeavor aiming for high accuracy, but 
correcting for bias is a also an important step in the construction of optimal parameter estimates. 
In this article, it is demonstrated how a popular model of confinement can be corrected for finite 
sample bias in situations where the underlying data exhibits Brownian diffusion and observations 
are measured with non- negligible experimental noise (e.g., noise induced by finite photon counts). 
The work of Tang and Chen [J. Econometrics^ 149 (2009)] is extended to correct for bias in the 
estimated "corral radius" (a parameter commonly used to quantify confinement in SPT studies) in 
the presence of measurement noise. It is shown that the approach presented is capable of reliably 
extracting the corral radius using only hundreds of discretely sampled observations in situations 
where other methods (including MSD and Bayesian techniques) would encounter serious difficulties. 
The ability to accurately statistically characterize transient confinement suggests new techniques 
for quantifying confined and/or hop diffusion in complex environments. 

PACS numbers: 87.80.Nj, 87.10.Mn, 05.40Jc, 2.50.Tt, 5.45.Tp 
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I. INTRODUCTION 

In many live cell applications, large scale cellular 
structures impose complex constraints on the motion of 
smaller biomolecules p^HT3]. Quantifying these effects 
from in vivo observations is the goal of numerous exper- 
iments. Fortunately, recent advances in microscopy and 
other single-molecule probes have substantially improved 
resolution in both time and space, so various complex ki- 
netic constraints can be more quantitatively measured. 

Fluorescence microscopy can be used to extract ki- 
netic information from a sequence of point spread func- 
tion (PSF) measurements [Mj |T5]. Pioneering efforts 
[I6l [17] aiming to quantify transient "corralling" param- 
eters characterizing confinement induced by cytoskeletal 
and transmembrane protein structures [3 appealed to 
Mean Square Displacement (MSD) analyses. Recently, 
the utility of more statistically motivated time series 
analysis have become more popular for analyzing single- 
molecule data. These tools offer several advantages over 
traditional MSD-based techniques. For example, like- 
lihood and Bayesian-based statistical analysis methods 
permit one more flexibility in terms of inference decisions 
characterizing noisy systems and these schemes also pro- 
vide more efficient estimation strategies [Sj [131 [IBHIS] . 

Many of the first works utilizing likelihood-based anal- 
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ysis methods to analyze single particle tracking (SPT) 
data ignored the effects of measurement noise, but 
the importance of modeling this noise source has been 
demonstrated in various works focused on analysis single- 
molecule data where measurement noise induced by the 
experimental apparatus is not negligible relative to ther- 
mal ffuctuations inherent to single-molecule measure- 
ments [22-25 . Ref. [23 provides a discussion on issues 
associated with simultaneously quantifying measurement 
and molecular diffusion in SPT applications, but the fo- 
cus of Ref. [23 is on optimal parameter estimation. It 
is well-known that the maximum likelihood estimator is 
asymptotically unbiased [18 and achieves the Cramer- 
Rao lower bound when the assumed underlying model 
precisely matches the data generating mechanism pro- 
ducing observations [T9]. In applications where tracking 
molecules for a long time is complicated due to crowd- 
ing, photobleaching, and/or emitter "blinking" [14", =15] , 
it is difficult to collect a large number of measurements 
(hence the asymptotic sampling regime is not encoun- 
tered). In PSF modeling, the appropriate parametric 
models have been more broadly agreed upon [19], but 
the "correct" stochastic model consistent with experi- 
mental single-molecule observations is a more delicate 
issue [25l[26]. Furthermore, even if observations are con- 
sistent with the assumed stochastic model, correcting for 
systematic bias introduced by finite sample sizes where 
observations contain both diffusive noise and measure- 
ment noise has not received great attention in the SPT 
literature. So estimators accurately quantifying finite 



sample bias (as opposed to asymptotically minimizing 
parameter variance or bias) are desirable when analyzing 
experimental trajectories. 

This work introduces a bias correction scheme for ex- 
tracting the "corral radius" [IGf jlT^. This quantity is 
commonly used to characterize confinement in biophys- 
ical applications [71 |Tl] |T3] . Examples characteristic of 
sampling regimes encountered in fluorescence microscopy 
are presented, but the approach can be readily general- 
ized to other time and length scales. The bias correction 
removes systematic errors induced by observing a short 
finite time series (this enables estimation in situations 
where an MSD curve is deemed too noisy for information 
extraction) in contrast to removing artifacts of motion 
blur |5] ^22^. However the analysis presented explicitly 
shows how to map estimated parameters to MSD curves, 
so previously proposed motion blur corrections for con- 
finement O |22] can be used to augment the parameter 
bias results presented. 

Likelihood-based techniques are employed throughout 
[2T| \25\ [27H3T] ; such methods enable one to consider nu- 
merous time series analysis tools in physical and life sci- 
ence applications. The author has found adopting sta- 
tistically rigorous time series analysis tools from econo- 
metrics and computational finance helpful in statisti- 
cally analyzing data from microscopic simulations [32- 
[34] and single-molecule force manipulation experimen- 
tal data where measurement noise is commensurate with 
thermal noise [24l [25l [35] . In this article the relevance of 
recent likelihood-based tools to SPT modeling is demon- 
strated. Section [ll] presents the stochastic differential 
equation (SDE) model considered, relates parameters ex- 
tracted from these models to traditional MSD analyses, 
and introduces the basic tools utilized throughout. The 



first figure and tables in Sec. Ill present the main re- 
sults; the remaining results explain and justify how a 
theory originally developed for monitoring SDKs without 
measurement noise [29 can be modified and extended to 
handle the situation where measurement noise contami- 
nates observations. 



II. METHODS 

The underlying position of a molecule will be denoted 
by X and the noisy experimental observations will be de- 
noted by y; the motion models considered take the fol- 
lowing form: 



dxt = - W{xt)dt + crdBt 



(1) 

:A/-(0,i?), (2) 



The above is an SDE model [36] with a constant diffusion 
coefficient driven by a standard Brownian motion process 
Bt (the subscripts denote a continuous time model) hav- 
ing a drift function determined by a potential V{x). The 
measurements i/i in Eqn. [2] are contaminated by noise e 
assumed to be random variables drawn from a Normal 



distribution with mean zero and variance R (in SPT, one 
often focuses on the effective measurement noise quan- 
tified by R^^'^); the measurement noise is assumed to 
be an independent and identically distributed (i.i.d) ran- 
dom number sequence and the variance R is assumed 
unknown a priori (the model also assumes statistical in- 
dependence of Xi and e^). The integer subscript i de- 
notes that trajectory observations are made at discrete 
times and t^+i — ti = At\/ i. Since typical SPT calcula- 
tions assume independence between spatial coordinates 
[SI [TTJ [171 123] , we will restrict attention to analyzing the 

ID version of Eqn. [ll hence the diffusion coefficient is 

2 — 

D •= — 

For V{x\ two different functional forms will be con- 
sidered: (i) V{x) = for \x\ < L/2 and V{x) = 00 for 
\x\ > L/2 which we refer to as reflected Brownian mo- 
tion (RBM); in this case the parameters needed to com- 
pletely characterize particle motion are (L, a, R) and (ii) 
V{x) = ^tvx'^ which we label as the Ornstein-Uhlenbeck 
(OU) process (also known as the Vasicek model [29]): 
parameters requiring estimation in this case are (/>:, a, R). 
Both potentials mentioned above have been considered in 
confined membrane diffusion studies [5l [131 III] • Kusumi 
et al. demonstrated how the MSD asymptotically ap- 
proaches ^ in the RBM model; extraction of L from 
data is still a common technique for quantifying confine- 
ment in SPT studies [5] [Til [13] (the ID corral radius is 

defined by \r^)- 

The MSD corresponding to an ergodic OU process ob- 
served with infinite time is (see Appendix): 



MSD((5) = — (1-e" 



-k5\ 



(3) 



If the two models under consideration have identical 
diffusion coefficients, then setting k = ^ is one way 
to match asymptotic MSD parameters; in the confined 
regime, this relation also allows one to map k, of the OU 
onto the corresponding L in the RBM model. The Ap- 
pendix displays representative trajectories and also com- 
pares the entire MSD for OU and RBM models driven 
by the same Brownian noise realization. In the measure- 
ment noise free case {R = 0) with large samples, the 
RBM and OU processes are easy to qualitatively and 
quantitatively distinguish. When measurement noise is 
present {R > 0), the two scenarios are much harder to 
distinguish if one only has access to a few hundred ob- 
servations of each trajectory. It is shown that for 100- 
400 observations, that hypothesis testing tools [3T can- 
not statistically distinguish the two models in parameter 
regimes of relevance to many SPT studies. If statistical 
signature of other more complex noise cannot be system- 
atically detected in these sample sizes [4, 7, 9, 12, 38 , one 
should consider modeling with the OU process because of 
statistical advantages this process offers when analyzing 
experimental data (these are discussed in the next sub- 
section). The advantages (from a physical standpoint) of 
applying detailed time series analysis to short trajecto- 



ries experiencing transient confinement are discussed in 
the Conclusions. 



A. Advantages Afforded by the OU Model 



The discrete time analog of Eqn. [T] for the OU model 



Xi =Fxi-i 
yi =Xi + e 



2KAt\ 



(4) 
(5) 



where F = e'^^^ and Q = |^(1 - e'^'^^^) [29 . This 
relation allows one to readily use the Kalman filter es- 
timation framework [27 . Maximum likelihood estima- 
tion (MLE) of the parameters completely characterizing 
the stationary OU process can be computed from the 
observable measurements {yi}fLo [24l|27l[35]. This per- 
mits efficient estimation in situations where sample sizes 
for an MSD analysis are difficult to reliably extract and 
statistically characterize (see Appendix Fig. IgI). The 
Gaussian structure of the OU process also enables one to 
exploit a variety of other powerful tools that can be used 
to analyze this type of stochastic process [27], includ- 
ing goodness-of-fit testing (checking model assumptions 
against data directly ^24r.26] ) . exact rate of convergence 
analysis under stationary and non-stationary sampling 
[28 , and bias correction. For example, Tang and Chen 
[29 demonstrated how to remove bias from MLEs com- 
puted using finite sample sizes in the case where the x^'s 
are directly observed (i.,e., R = 0). In the stationary 
case {hz > 0), it can be shown using moment bounds for 
weakly dependent sequences (more specifically a— mixing 
random variables [29l |39l [40] ) that: 



E[k] 



(6) 



1 






^At 
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Oij,), 



NAt^2 ' ^ '2^ 

where E[k] denotes the expectation of the estimate k. 
The other terms quantify the expected bias induced by 
finite N (time series sample size) . k is often the dominant 

source of bias when the relation L = \/^^— is used to 
extract the corral radius from OU parameter estimates 
in the sampling regimes studied (e.g., results obtained by 
plugging in the corrections to <j^ reported in Ref. \2W did 
not affect results). 

Before moving onto the case where i? > 0, it is worth 
reviewing a classic first order autoregressive time series 
model [27, 28 where Xi = Fxi-i + t] where 77 ~ A/'(0, Q); 
the interest is in estimating F (Q is considered frozen and 
to be nuisance parameter). For notational simplicity set 
(5 = 1 and xo = 0. In this case, the standard likelihood 
equation is: 

1 ^ 
Pf{xi,X2,...,xn) = (27r)"^exp(- ^^{^i - i^^i-i)^) 

i=l 

(7) 



Taking the logarithm, expanding the quadratic terms, 
and setting the derivative of the expression above with 
respect to F equal to zero provides the classic MLE esti- 
mate: 



F 



N 



(8) 



The above suggests a naive suboptimal (denoted by a 
tilde) estimator: 



N 

E(^i + ei)(xi_i +ei_i) 

i=l 

N 



(9) 



where R is an independent estimate of the measurement 
noise variance. Recall that the measurement noise is as- 
sumed i.i.d., so if R is asymptotically consistent and Q is 
fixed, F is asymptotically consistent since the cross-term 
sums involving e and x tend to zero and become insignif- 
icant relative to the other non-zero sum in the k, > case 
under study. The problem with this approach is that the 
estimator is suboptimal (the cross-terms increase estima- 
tion variance). Unfortunately, the estimator above also 
requires one to construct a consistent R (this can alter- 
natively come from a prior, but this will likely introduce 
bias which is hard to quantify). Furthermore, if one uses 
estimators ignoring confinement effects, new systematic 
biases (on top of inherent finite sample bias associated 
with estimating Hi) can be introduced. This phenomenon 
is demonstrated by example in the Results. 

In the Kalman filter and ID OU model considered, 
the innovation likelihood (Appendix Eqn. [TT] ) has an ap- 
proximate autoregressive [27 form if the filter covariance 
reaches steady state quickly (this occurs in the param- 
eter regimes explored). If a stationary OU process is 
deemed adequate to describe experimental observations 
and the Kalman filter covariance sequences reaches its 
steady state value rapidly, then analysis in Ref. [29] can 
be applied to study the expected finite sample bias of k. 
When one jointly estimates the MLE parameters associ- 
ated with Eqn. [ijby optimizing the innovation likelihood 
(see Appendix Eqn. [TT] ), one effectively returns to the 
situation in Eqn. [S] where the estimates of F can be ex- 
tracted without knowledge of the value of the constant 
noise parameters. Note that when the OU process is dis- 
cretely observed without noise {R = 0) and parameters 
are estimated jointly, it can be proven that tz and a^ the 
MLEs of k and a^ are asymptotically independent 09] 
(i.e., estimation results of k are not affected by a^). In 
the Results (Fig. [4|, the convergence of matrices char- 
acterizing the Kalman filter are demonstrated (in con- 
fined diffusion, these matrices typically converge rapidly 
to steady-state). 

In what follows, it is shown how plugging the MLE's 
Eqn jT into Eqn. [6] can significantly reduce bias from 



parameter estimates associated small N in parameter 
regimes of relevance to SPT tracking (the approach 
avoids specifying the "lag parameter" plaguing MSD- 
based analyses [2T). The approach is demonstrated to 
accurately infer both Hi and the effective L (corral ra- 
dius) if data is generated using either the OU model 
(correct model specification) or the RBM (model mis- 
specification). 



III. RESULTS 

Figure [T] presents a histogram of the raw estimate of 
the corral radius obtained via the relation L = \/^^ 
for the case where 1000 Monte Carlo simulations with 
L = 400nm,D = 0.2/j.m'^ /2,R^^'^ = 50nm, At = 25ms 
and A^ = 100 observations of y were used to generate 
data. From this data, the MLE parameter estimates 
characterizing the model in Eqn. [ijwere extracted. Cor- 
ral radius parameters are inferred using the MLE and the 
bias corrected parameter estimates for two different data 
generating processes. In the top panel, the OU process 
generates data; in the bottom panel, the RBM process 
generates data (here there is model mismatch). The bias 
induced by only observing 100 time series is effectively 
removed in both cases. Appendix Fig. [6] displays rep- 
resentative trajectories of x, y, and the empirical MSB 
associated with these trajectories. 
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TABLE I. Corral radius estimates with innovation MLE and 
bias corrected MLE. The columns labelled with L contain the 
average parameter estimate obtained by analyzing 200 Monte 
Carlo trajectories each containing 400 observations spaced by 
At = 25 ms (the number in parenthesis reports standard 
deviation) . The column labeled error reports the mean minus 
known true corral radius. In this table D = 0.2/im^/s and 
i?i/2 = 25nm. 





OU 


RBM 


Estimator 


L [nm] 


Error 


L [nm] 


Error 




L = 250nm 


Classic Innov. 


169.31 (16.07) 


-80.69 


168.19 (12.56) 


-81.81 


Bias Cor. 


241.40 (22.54 ) 


-8.60 


240.24 (17.56 ) 


-9.76 




L = 400nm 


Classic Innov. 


277.09 (18.38) 


-122.91 


275.69 (10.63) 


-124.31 


Bias Cor. 


399.09 (26.89 ) 


-0.91 


398.21 (15.39 ) 


-1.79 




L — 500nm 


Classic Innov. 


344.25 (25.75) 


-155.75 


343.70 (14.83) 


-156.30 


Bias Cor. 


500.35 (38.69 ) 


0.35 


501.73 (22.54 ) 


1.73 



Tables |ipi| present similar results, but vary the system 
and sampling parameters. The parameters explored were 
motivated by SPT studies. Even for N = 400, substantial 
bias exists in the asymptotically efficient MLE. Bias re- 
duction comes at the cost of variation as can be observed 



FIG. 1. Raw MLE and bias corrected corral radius estimate, L, 
obtained by extracting the parameters from time series of length 
100 sampled every 25ms (this was repeated for 1000 Monte Carlo 
trials; the histogram displays 1000 L's). In the top panel, the OU 
process (with known parameters) generates observations. In the 
bottom panel, the same Brownian motion paths used to generate 
OU trajectories are used to construct RFM paths. In both cases, 
a single measurement noise random number stream was added to 
each trajectory. The use of the same underlying Brownian path and 
measurement noise sequence was used to reduce random variation 
and facilitate quantifying systematic errors. 



by the reported standard deviations. However, using the 
OU model structures allows one to use a wealth of quan- 
titative tools for understanding experimental data anal- 
ysis. The main results have now been presented, what 
follows expands on technical details and on the domain 
of applicability of the bias removal approach. 

Figure [2] presents the distribution of k for a variety 
of estimators. The focus is on hi since this is often the 
major source of variation in L estimated from short time 
series [29 . The top panel displays three estimators; (i) 
the MLE associated with the OU process where R is as- 
sumed zero [29]; (ii) the MLE obtained by jointly ex- 
tracting the hv^a'^, and R that minimize the innovation 



likelihood [35] (see Eqn. 11) and; (iii) using the bias cor- 
rection of Tang and Chen [29 applied to the output of 
(ii). The average of the k distribution for the three cases 



TABLE II. Same as Table [l] except D = 0.02fim'^/s. 





OU 


RBM 


Estimator 


L [nm] 


Error 


L [nm] 


Error 




L = 250nm 


Classic Innov. 


169.43 (19.17) 


-80.57 


170.70 (12.24) 


-79.30 


Bias Cor. 


253.54 (31.84) 


3.54 


258.49 (21.06) 


8.49 




L = 400nm 


Classic Innov. 


253.84 (46.88) 


-146.16 


257.78 (34.79) 


-142.22 


Bias Cor. 


418.74 (118.61) 


18.74 


437.48 (90.10) 


37.48 




L = 500nm 


Classic Innov. 


310.85 (61.59) 


-189.15 


308.72 (57.53) 


-191.28 


Bias Cor. 


582.40 (220.27) 


82.40 


603.12 (244.42) 


103.12 



are 24.3, 18.2, and 16.0 ^, respectively (the true value is 
15 -). The difference may seem small, but recall that the 
estimated corral radius depends nonlinearly on k (hence 
the amplified difference in L). 

The bottom panel in Fig. [2] uses "other" suboptimal 
estimators of Hi. In one case R is assumed to be known 
accurately a priori; in this case we set R^^'^ = O.SR^^'^ = 
AOnm and used Eqn. [9] to estimate the parameter. Since 
accurate a priori knowledge of R is can be a question- 
able assumption in SPT studies, we also show results of 
applying Eqn. [9] in conjunction with the estimator re- 
ported in Ref. ^22] to extract R from the data; here R 
is biased because the estimator in Ref. [22] was designed 
for the case where no forces or confinement constraints 
affect particle dynamics (the average of the MLEs assum- 
ing the model in Ref. [22] was 72.4 nm for this data set; 
note that Refs. [22", '2T warn that the estimator is not 
valid if constraint forces are present). 

Note that the average MLE (without bias correction) 
for the diffusion and measurement noise was (0.23 (im? / s^ 
41.9 nm), that for Ref. [22 was (0.07 iirn^/s, 71.3 nm), 
and the true value for the OU data generating process 
was (0.20 jim? / s, 50.0 nm). The arguments used to ex- 
plain the validity of the bias correction of Ref. [29] in 
conjunction with the Kalman filter's innovation likeli- 
hood [25^, 1271 (relevant expressions shown in Eqn. 11) 



are not directly applicable to the bias correction of the 
other parameters reported in [29^ when i? > 0. Analysis 
of the bias and variance of parameters a and R are more 
involved due iteration introduced by the Kalman filter 
update and forecast steps and this analysis is beyond the 
scope of this work. 

Application of various estimators of OU parameters to 
data generated by both the OU and RBM (a misspecified 
model) processes, was carried out for two reasons: (i) to 
emphasize that certain estimators (or Bayesian priors) 
can induce subtle systematic biases and (ii) stress that 
likelihood-based inference permits other analysis tools 
beyond estimation. Bias correction is possible in addi- 
tion to other techniques. For example, detecting con- 



finement from observations using visual inspection of the 
short trajectories is problematic (Fig. [6] shows how even 
in the i? = case, distinguishing RBM from the OU pro- 
cess is difficult with N = 100) . However, goodness-of-fit 
testing can be employed [26l ^\. Applying the tech- 
nique of Hong and Li [37] (more specifically computing 
the M(l, 1) test statistic) allows one to reject ~ 20% of 
the trajectories assuming the so-called directed diffusion 
model (i.e., constant diffusion, measurement noise and 
velocity [TT.'ir, but hz = 0) even with TV = 100. There is 
overwhelming statistical evidence for all TV = 400 cases 
(the average p— value obtained assuming the directed dif- 
fusion plus measurement noise model was < 5 x 10~^ for 
all N = 400 cases considered). This demonstrates that 
the test has power to detect kinetic signatures of con- 
finement in the presence of diffusive plus measurement 
noise in regimes of interest to SPT studies (if the model 
wasn't rejected one can safely use other estimators e.g., 
[22| |23]). The case where we assumed an OU model, 
but an RBM model actually generated the data (model 
misspecification), was statistically indistinguishable us- 
ing tests in Ref. [37] from the case where the OU model 
generated data. Since there is no evidence in the raw ob- 
servational data favoring one model over the other, and 
both models produce similar estimates of the quantity 
of interest (the corral radius), it is attractive to use the 
OU modeling viewpoint since a substantial body of liter- 
ature exists for analyzing this type of stochastic system 

[miiiiiiTHai]. 



Figure |3] provides another example of analysis tools 
that are made available from likelihood-based analyses. 
Here hz is plotted against the truncated bias expansions 
in Eqn. [g] (taken from Tang and Chen [29 ) for a fixed At 
and two sample sizes N. Note how as hi decreases, the 
fraction of bias increases rapidly. Also note, that as hi 
decreases, there is a higher likelihood of an MLE param- 
eter estimate being near or less than zero (even for a truly 
stationary process where the underlying data has z^: > 0) . 
A high value of hi suggest weak "corralling" since L is in- 
versely related to k. The inverse dependence also causes 
the inflated standard deviation for L = bOOnm since a 
small fraction of estimated k are near zero (also note that 
the median L's corresponding to the L = 500nm row of 
Tab. [TTJwere 540.3 and 582.4; this suggest that these esti- 
mates in the tail of the estimated parameter distribution 
substantially influenced the observed mean). Although 
one can remove expected bias, if the fraction of bias is 
large relative to the signal other factors can complicate 
bias correction. For example, (i) higher order terms in 
the expected bias expansion can become more important; 
(ii) inherent parameter uncertainty in the point estimate 
substantially affects the expected bias. Therefor plots 
like Fig. [3] allow researchers to to quantitatively deter- 
mine when other factors influencing the bias correction 
scheme need to be considered. 
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FIG. 2. Distribution of k obtained when the OU process with 
measurement noise generated observation (results correspond to L 
histogram shown in Fig. fTJ. Two different classes of estimators 
were used: (i) MLE-based estimators use a likelihood with the 
correct model and (ii) "Other" estimators use an assumed prior 
input for R (the true R^^"^ is 50 nm, but R}^"^ is set to 40 nm since 
knowledge of the precise effective measurement noise is difficult 
to accurately quantify 22 in SPT applications) and a suboptimal 
estimator in Eqn. l9l (for R we plug-in the noise estimate obtained 
using code associated with Ref. [23]). In both cases, the Innovation 
MLE (without bias correction) serves as the reference histogram. 



A. Conditions Required for Bias Removal 

The Kalman filter's constant noise assumption (i.e., 
the covariance of the innovation sequence, 5^, in Eqn. 
11) needs to be tested in order for the analysis of Tang 
and Chen [29] to be accurate for the expected bias in k. 
Figure [4] illustrates that this is indeed the case for the pa- 
rameter regimes under study. Note that we intentionally 
ignored the estimation of the mean of the OU process 
(the mean was set to zero), this simplifies analyzing the 



effect of n on the autoregressive parameter (F = e '^ ) 
under the assumption of a constant innovation covari- 
ance. The mean zero OU process does not restrict utility 
(with extra work one can analyze the joint mean and n 
estimates and in practice one can simply subtract the 
empirical mean of y). 
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FIG. 3. Bias vs. tz for two different sample sizes A^ using Eqn. [6] 
(At = 2hms and measurement noise magnitude 25nm). 



Beyond testing constancy of the covariance of the inno- 
vation, one should verify that the bias in the estimated R 
and cr'^ is small in relation to the bias in the estimated K:. 
This can be achieved via simulation if need be. If bias in 
R and a^ is determined to be significant, new analytical 
or numerical bias removal schemes should be considered 
(e.g., parametric bootstrap methods [29 [ [42 j ). 

When At is decreased with R fixed, measurement 
noise often becomes a more dominant part of the single- 
molecule signal and the bias in R is relatively small 
[m [30]. In the OU model, the influence of a^: on Q is 
0(At2) since Q = ^{2i^^t + 0{At^). In the small At 
limit, one can leverage existing nonparametric tools for 
estimation and inference, e.g. [31^. However, the param- 
eter regime explored here is one in which /^'s influence 
is not small relative to At (otherwise, the approach in 
Ref. [23 would predict more accurate R estimates). In 
this study, it was empirically demonstrated that the bias 
connected to the innovation MLE of R is small relative 
to that of K, in several parameter regimes of relevance to 
SPT modeling (no bias correction was applied to R). 

Note that the bias correction scheme presented also 
depends heavily on the stationarity assumption of both 
the state and on the innovation sequence. If stationarity 
of X is questionable, computing a corral radius should 
be reconsidered. More formally, unit root tests [27] (ac- 
counting for measurement noise) can be used to check the 
Brownian motion vs. stationary OU models. To more 
generally test stationarity of the mean or covariance of 
the observed measurements, other testing procedures can 
also be considered, e.g. [43 . Even if all stationarity tests 
are passed, if the (bias corrected) estimate and the as- 
sociated parameter uncertainty suggest k is near zero, 
the suitability of a confined diffusion model needs to be 
carefully reevaluated. 

Finally, if all conditions mentioned in this subsection 
are met and the Kalman filter corresponding to the OU 
process is an adequate model of the observations (an as- 
sumption tested here with time series hypothesis testing 



[37]), then the state and innovation noise residuals have 
mean zero (these residuals makeup stationary process un- 
der the conditions above). The filter and measurement 
noise sources make the innovation sequence different than 
classic order one autoregressive process, but the addi- 
tional noise terms do not substantially influence the first 
order expansion of the expected bias of k. The effects of 
the additional noise terms are are lowest when the scalar 
(see Eqn. 11), is close to one (a standard au- 



gam, 



K, 



toregressive process generates the data when Ki = 1). In 
small time series sample sizes, even when i^^ < ^, the 
bias correction is accurate since the effects of parameter 
uncertainty tend to dominate the additional noise asso- 
ciated with the filtered state estimates. Furthermore, 
the MLE parameters are found by jointly optimizing the 
Eqn. [TT] given data, but when the innovation covariance 
quickly reaches steady state, the analysis of Tang and 
Chen [29] is relevant. 
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FIG. 4. Innovation covariance (see Eqn. |ll| ) convergence rates. 
Plots were obtained by plugging in exact data generating parame- 
ters studied in Tables |l|lT] The plot illustrates that the filter quickly 
reaches steady state (total sample sizes ranged from 100-400 obser- 
vations). 



IV. CONCLUSIONS 



sizes considered, there was not adequate evidence to dis- 
tinguish reflected Brownian motion from an OU process 
observations alone. The estimated corral radius was re- 
liably extracted using an OU model regardless of the 
underlying stochastic dynamics and measurements pro- 
ducing the data. Numerous statistically motivated rea- 
sons for favoring the OU model to the reflected Brow- 
nian motion model were discussed. Potential problems 
that can be encountered when the inferred corral radius 
is too large to reliably infer from the data available were 
also discussed. The rich likelihood structure afforded by 
wrapping SDE plus noise models around experimental 
data was exploited throughout. The likelihood formula- 
tion circumvents the the need for selecting ad hoc sam- 
pling parameters such as a "time-lag" cut-off (this is a 
common problem MSD-based analysis [23 which are still 
quite common and popular in the SPT community). 

The ability to accurately extract kinetic parameters 
and correct for biases induced by small time series sample 
sizes (while also accounting for measurement and thermal 
noise in a statistically rigorous fashion [2?, ^25^) shows 
great promise studies where the underlying molecule ex- 
periences random forces whose distribution changes in 
both time and space due to complex interactions in a 
highly heterogeneous environment. For example, if one 
can both reliably determine when molecules leaves a 
"picket fence" [3 in the plasma membrane via change 
point detection algorithms [44 and can track trajecto- 
ries with high temporal resolution (perhaps at the cost 
of spatial accuracy), one can utilize the tools presented 
here to accurately map out both the diffusion coefficient 
and the corral radii explored by molecules in the plasma 
membrane. This presents an attractive physically inter- 
pretable modeling alternative to sub-diffusion or contin- 
uous time random walk type models, but such a study is 
left to future work. The method introduced was shown 
to be useful in parameter regimes commonly encountered 
in fluorescence-based SPT experimental studies, but the 
approach is general and can be used to probe other length 
and time scales afforded by the measurement device. 



Simulations were used to demonstrate that kinetic con- 
finement parameters could be accurately extracted from 
relatively small sample size (100 < N < 400) time se- 
ries containing both inherent diffusive and measurement 
noise (the latter prevents direct observation of position) . 
A bias correction scheme expanding off of Ref. [29] was 
presented; it was demonstrated that the scheme can ac- 
curately extract the corral radius in parameter regimes 
commonly encountered in membrane diffusion studies. 
The reasons why the approach work were discussed. 

Two popular data generating processes were consid- 
ered (reflected Brownian motion and the OU process). It 
was demonstrated that accurate results can be obtained 
even if the stochastic model assumed was not consistent 
with the data generating process providing some robust- 
ness assurance. In the confinement regime and sample 
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VI. APPENDIX 



A. MSD of the Stationary {n > 0) OU Process 



The MSD of a generic process is MSD ((5) := 

T 

^(Z](^t+5 — Xt)'^). Here (•) denotes ensemble averaging 
[SI HI] ; plugging in the solution to the mean zero station- 



ary OU process (variance = |^ [36^) and exploiting other 
standard properties of SDEs driven by Brownian motion 
[45] , one can derive the fohowing: 
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FIG. 5. Long time sample trajectory for OU and RSB process 
with (top left) and without (top right) measurement noise. In this 
case, the sample matches the theoretical MSD limit nicely. Here 
L = 400nm, D = ^.2^m? /2, R = 50nm, and At = 25ms. 

To account for i.i.d. Gaussian measurement noise (i.e. 
one carries out an MSD on y) in the above expression, 
simply add T2R to the MSD expression above |3Qj . 



B. Representative Trajectories and MSDs 

In this section, the reflected Brownian motion and the 
corresponding OU process (found using Eqn. [6| are plot- 
ted with and without measurement noise. The MSDs of 
the measurement noise free and measurement noise case 
are shown for both large and small sample sizes. 
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FIG. 6. Same as [s] except the sample size was reduced to 100 
observations. The small sample sizes complicates reliably using 
an MSD-based analysis. In small sample sizes, well-known issues 
associated what lag truncation to use, statistical dependence com- 
monly introduced when computing MSDs, etc. |23| are even more 
pronounced. The likelihood-based bias correction scheme intro- 
duced is able to reliably extract system parameter even with these 
small samples sizes. 



C. MLE of the Innovation Sequence 

In the main text, mappings between the OU param- 
eters and those of the classic Kalman filter [2T were 
presented. Here the equations defining the innovation 
MLE and the associated likelihood [24l [27] relevant to 
the scenario studied are presented (the reader is referred 
to the excellent monograph [27 for details). Note that 
the state here is x and the "observation matrix" H and 
that yi\i-i = Hxi\i_i (= Xi\i_i in the case considered). 

(11) 
{R, F, Q) = argmax C{R, F, Q) = p{yi,y2, ...,yN]R, F, Q) 

p{yi,y2,'",yN]R,F,Q) = 



N 



n — 172 ^^p(" 



-JVi -Fyi-i\i-iy 

2S, 



yi\i = yi\i-i + Ki{yi - yi\i-i) 
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For the stationary OU process, the recursion above 
(processing the observation) was started using f i|o = 



and Pi|o = f^. The Nelder-Mead algorithm was used to 



find the parameter optimizing Eqn. 11 Goodness-of-fit 



testing [261 El] was used to both check the consistency 
of model assumptions against data and to ensure that a 
local minimum was not encountered in the optimization. 
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